Numerical simulations of strong incompressible magnetohydrodynamic 
turbulence 

J. Mason, ^ J.C. Perez, ^ S. Boldyrev,^ and F. Cattaneo^ 

Department of Astronomy & Astrophysics, University of Chicago. 5640 S. Ellis Ave, Chicago, IL, 60637, 
USA 

Institute for the Study of Earth, Oceans, and Space, University of New Hampshire, Morse Hall, 8 College Road, 
Durham, NH, 03824 

Department of Physics. University of Wisconsin at Madison, 1150 University Ave, Madison, WI 53706, 
USA 

(Dated: 17 February 2012) 

Magnetised plasma turbulence pervades the universe and is likely to play an important role in a variety 
of astrophysical settings. Magnetohydrodynamics (MHD) provides the simplest theoretical framework in 
which phenomenological models for the turbulent dynamics can be built. Numerical simulations of MHD 
turbulence arc widely used to guide and test the theoretical predictions; however, simulating MHD turbulence 
and accurately measuring its scaling properties is far from straightforward. Computational power limits the 
calculations to moderate Reynolds numbers and often simplifying assumptions are made in order that a wider 
range of scales can be accessed. After describing the theoretical predictions and the numerical approaches that 
are often employed in studying strong incompressible MHD turbulence, we present the findings of a series of 
high-resolution direct numerical simulations. We discuss the effects that insufficiencies in the computational 
approach can have on the solution and its physical interpretation. 
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I. INTRODUCTION 

Magnetohydrodynamic turbulence constitutes one of 
the most important unresolved problems in classical 
physics. Magnetohydrodynamics provides the simplest 
theoretical framework in which we can develop our un- 
derstanding of magnetised plasma turbulence, which it- 
self forms the foundation for a vast array of astrophysical 
phenomena that are believed to be magnetically driven. 
However, despite it being almost 50 years since the pio- 
neering work of Iroshnikovi & Kraichnan^, debates con- 
tinue over even the most fundamental issues, such as the 
inertial range scaling properties of the energy cascade. 

We consider here the simplest case of incompressible 
field-guided MHD turbulence. It is useful to write the 
equations governing the evolution of the fluctuating ve- 
locity v(x, and magnetic field b(x, in terms of the 
Elsasser variables = v ± b, 

(^^ T Va • z± + (zT • V) z± = -VP + z/V^z^+f^, 

(1) 

V • z± = (2) 

where Va = "Bq/ ^/4^Tpo is the Alfven velocity based 
on the uniform background magnetic field Bq, pq is the 
background plasma density that we assume to be con- 
stant, P = {p/po + b'^/'2) is the total pressure, i/ is the 
fluid viscosity (which for simplicity has been taken to 
be equal to the magnetic diffusivity) and f ^ represents 
random forces that drive the turbulence at large scales. 
The pressure term can be eliminated in favour of a pro- 
jection of the solution onto the 'incompressible plane' 



(see, e.g., Ref. Q). In the absence of forcing and dissipa- 
tion the linearised system admits solutions in the form 
of Alfven waves that travel parallel or antiparallel to the 
background field Bq = B^e^, say, with the Alfven speed. 
A normal mode analysis reveals the dispersion relation 
a;^(k) = ±fc|| V^. The waves are transverse and they can 
be divided into two classes: shear Alfven waves with po- 
larizations perpendicular to both Bq and to the wavcvec- 
tor k, and pscudo Alfven waves with polarizations in the 
plane of Bq and k, perpendicular to k. 

Interactions between counter-propagating Alfven wave 
packets transfer energy to smaller scales (Ref. [2j) until 
eventually the dissipative scales are reached and energy 
is removed from the system. The efficiency of the trans- 
fer is controlled by the relative size of the linear and 
nonlinear terms in equation ([T}. The regime in which 
the linear terms dominate is known as weak MHD tur- 
bulence, otherwise the turbulence is called strong. In 
fact, it has been demonstrated both analytically and nu- 
merically that the energy cascade occurs predominantly 
in the plane perpendicular to the guiding magnetic field 
(Refs.|3and[5|), which ensures that even if the turbulence 
is weak at large scales it encounters the strong regime as 
the cascade proceeds to smaller scales. We concentrate 
here on the strong case. 

The first phenomenological theories for the energy cas- 
cade in strong incompressible MHD turbulence were de- 
voted to the so-called balanced case, in which the energies 
in each of the Elsasser fields = i / {z^Y<Px are ap- 
proximately equal and hence there is no net cross helicity 
in the system, H" = /(v • b)^^^ = E+ - E- 0. For 
this case, Goldreich & Sridhar— argued that as the cas- 
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cade proceeds to smaller scales the linear and nonlinear 
terms establish a so-called 'critical balance' 

VAkii ^ zfk±. (3) 

Consequently, the nonlinear interaction time t ^ X±/z^ 
balances the Alfven time ta ~ \\/Va and it follows that 
the total field-perpendicular energy spectrum E{kj_) oc 

— 5/3 

kj_ , with the field-parallel and field-perpendicular 

lengthscales related by A|| ex A^'''^ (Ref. 0). Although nu- 
merical simulations at the time verified the anisotropic 
cascade (e.g. Refs. d, 0-0), doubts began to surface 
in later years when higher resolution simulations with 
stronger guide fields seemed to be more consistent with 
E{kA_) oc kl^^^ (e.g. Ref. i, IMl). An explanation 
was provided by the theory of scale-dependent dynamic 
alignmenti^iii. The theory predicts that the fluctuating 
velocity and magnetic fields align within a small angle 
9 oc A^^^ in the plane perpendicular to the background 
field. Alignment reduces the size of the nonlinear term by 
an amount proportional to 9 and yields E{k±) oc kj^^^. 
Numerical simulations of driven, globally balanced, field- 
guided MHD turbulence have verified the alignment scal- 
ing^i^i^ and have shown that the domain fragments 
into regions of highly aligned and anti-aligned velocity 
and magnetic field fiuctuation a^^i^^ . Thus, even in the 
globally balanced case, cross-helicity plays a crucial role 
locally. 

Recently, attention has moved on to the globally un- 
balanced (or cross-helical) system in which ^ E^ . 
In this case the timescales for the nonlinear deforma- 
tion of the wavepackets ~ X±/z^ can be con- 
siderably different. There are a number of competing 
theoretical predictions that differ in regards to what 
is assumed about the physics of the nonlinear cascade. 
For example, the two different theories by Lithwick et 
alJ^ and Beresnyak & Lazarian^S ultimately arrive at 
the same conclusion that in the unbalanced regions the 
field-perpendicular Elsasser spectra have the same scal- 

ings E^{k±) oc kj_^^^. The two different theories by 
Perez & Boldyre\J^ and Podesta & Bhattacharjeo^i that 
are based on scale-dependent alignment propose that 
E^{k±) oc kj_^^^, while the theory by Chandran^^ con- 
cludes that the two Elsasser spectra have different scal- 
ings depending on the degree of imbalance. 

The goal of numerical simulations is to help clarify 
the picture. The aim is to measure the inertial range 
scaling properties and to use the results to discriminate 
between the conflicting theoretical predictions. However, 
although the task appears to be straightforward, the in- 
ference of scaling laws from numerical data is fraught 
with problems. The two main difficulties result from 
the difference between the theoretical predictions for 
the spectral exponents being very small, and from the 
fact that even state-of-the-art technological resources still 
limit the extent of the inertial range to approximately a 
decade. It is therefore of paramount importance to be 
able to make accurate numerical measurements, to en- 



sure that any source of error in the numerical data is 
minimised, and to invest all of the computational power 
in reaching an inertial range that is as extended as possi- 
ble. In the next section we describe a number of aspects 
of the simulation design and the techniques for data anal- 
ysis that enable us to measure the inertial range char- 
acteristics as accurately as possible with the currently 
available computational power. 



II. COMPUTATIONAL APPROACH 
A. Simplified equations 

Considerable simplifications can be made by making 
use of the structure of field-guided MHD turbulence. 
Since strong MHD turbulence is dominated by fluctua- 
tions with k± ^ , the polarization of the pseudo Alfven 
waves are closely aligned with Bq. Since field-parallel 
gradients are small, Goldreich & Sridhar— argued that 
the pseudo Alfven modes are likely to be dynamically 
insignificant. Indeed, if in equation ([T]) we neglect the 
term (z|| • V||)z^ in comparison with (z± ■ V^)z^ then 
the equations for the shear Alfven dynamics decouple 
from the pseudo Alfven dynamics and we obtain 

(§I^^A- V||^ z± + (z^ . Vi) z± = -V^P+uV'z^+i^, 

(4) 

V • zj = 0. (5) 

Equations (|4l5p are equivalent to the reduced MHD 
(RMHD) model that was originally derived in the context 
of fusion devices by Refs. [l^ and [13 (see also Ref. [25l ). 
We note that while this system has only two vector com- 
ponents = (z^,z^,0) each component is a function 
of all three spatial coordinates, x,y and z. Indeed, as 
stated in the critical balance condition the linear 
term involving field-parallel gradients balances the non- 
linear term involving field-perpendicular gradients. The 
RMHD system is therefore fundamentally different than 
the two-dimensional system in which dz = 0. 

The reduced MHD model has been widely used in 
the recent literature for studying the characteristics of 
strong field-guided MHD turbulence (see, e.g., Refs. [26l- 
[28[ ). Compared with system ([Ij, computing the Fourier 
series solution to system (|4]) is approximately twice as 
computationally efficient. Such savings can then be in- 
vested in reaching the higher Reynolds number regime. 
In mil Al below we directly compare the solutions of the 
RMHD model and the full MHD equations. We verify 
and put on firm ground that reduced MHD accurately 
captures the field-perpendicular dynamics of the strong 
turbulent cascade. 
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B. Parameter Regime 

We aim to simulate strongly nonlinear, field-guided 
MHD turbulence in which Vrms ~ ^rms Bq. If we 
choose the amplitude of the driving so that Vrms ~ 1: 
our previous work has shown that the universal proper- 
ties then set in when i?o ^ 3 (e-g- Ref- fl5l )- In all of the 
simulations reported in this paper we have set Bq = 5. 

Although weak turbulence eventually becomes strong 
as the cascade proceeds to smaller scales, failing to es- 
tablish the critical balance condition at the driving scales 
lengthens the transition region and therefore shortens the 
inertial range. An efficient way of driving strong turbu- 
lence while maintaining an inertial range that is as ex- 
tended as possible is to elongate the box in the direction 
of the guide field and to drive the lowest field-parallel 
and field-perpendicular wavenumbers fc|| = 27r/L||, k±^ = 
2tt/L±. Equation ([3|) is then satisfied at the forcing scales 
provided that 

^ zf/Bo (6) 

Here we fix L± — lix. In the balanced case in which 
the turbulence is driven so that z+ ~ ~ 1, the above 
condition can then be satisfied by taking L\\ ^ BqL±. 
In fact, the RMHD description already assumes that 
Bo ^ brms and therefore all RMHD simulations with 
L|| = BqLj^ produce equivalent results. 

The random forcing mechanisms constitute indepen- 
dent driving of both Elsasser populations, which allows 
us to vary the fluxes independently and study both 
the balanced and imbalanced regimes. As has been well 
documented in the literature, the particular mechanism 
of large scale driving is not essential for the scaling of 
the inertial interval (see, e.g., Refs. [lol and |l2h . Our 
driving mechanism for the MHD case has no component 
along z and is solcnoidal in the x — y plane. The ran- 
dom values of the Fourier coefficients of the forces in- 
side the range of wavenumbers 1 < |fcx|,|A;y| < 2 and 
(27r/L||) < \kz\ < n||(27r/L||) with n\\ = 1 or 2 are re- 
freshed independently on average every Q.\L^/ {2T:Vrms) 
time units (i.e. the force is updated approximately 10 
times per large-scale turnover time) with the amplitude 
of the force being chosen so that the resulting rms ve- 
locity fluctuations are of order unity. The variances 
^± = control the average rates of injection into 

the z+ and z~ fields. We take cr+ > (t_ and in the statis- 
tically state we measure the degree of imbalance through 
the parameter cr^ = / E = {E+ - E')/{E+ + £-). 
Thus balanced turbulence corresponds to CTc = 0, while 
the maximally aligned/imbalanccd case corresponds to 
(Tc = 1. 

The numerical resolution and Reynolds number are 
intricately related and one would like to conduct sim- 
ulations with as high a resolution and Re as possible. 
The upper bound on the resolution is determined by the 
availability of computing resources, with a doubling of 
the number of mesh points in each of the three direc- 
tions resulting in an increase in the computational cost 



by a factor of 16. For a chosen calculation size, the op- 
timal Reynolds number Re = Vrmsh/'^ ^ l/v can be 
established via a convergence study. Once a satisfactory 
Re has been established at low resolution, an estimate 
of the permitted value for higher resolution studies can 
be obtained from the scaling for the number of mesh 
points N ^ Re^ where /3 = 2/3 for the spectral expo- 
nent p = -3/2 (/3 = 3/4 for p = -5/3). A further note 
is required concerning the field-parallel resolution: Since 
the turbulent cascade is anisotropic we also allow for the 
numerical grid to be anisotropic, i.e. the spatial discreti- 
sation is performed on a grid with a resolution of N'j_ x A^|| 
mesh points. For cases in which A^|| is much less than N±_ 
we replace the diffusion operator in equations ([T]) and @ 
with v{dxx + dyy) + v\\dzz- 

We note that hypcrdiffusion is sometimes used to arti- 
ficially extend the inertial range, i.e. the diffusive terms 
in equations ([T]) and ^ are sometimes replaced with the 
higher-order operator (— l)"^^i/„V^" for n > 1. In hy- 
drodynamic turbulence, hypcrdiffusion is known to lead 
to a bottleneck effect (or a pile up of energy at large 
wavenumbers) that can ultimately affect the inference 
of scaling laws within the inertial range. Its effects on 
the inertial range dynamics for MHD are not well un- 
derstood, with simulations by different groups produc- 
ing different rcsultSiSSi^. Since our aim is to conduct as 
clean a simulation as possible, and since our task of dif- 
ferentiating between two very similar scaling exponents 
is difficult enough as it is, we prefer not to incorporate 
hypcrdiffusion into our simulations. 

Finally, it is important to have a large enough statisti- 
cal ensemble from which averages will be computed. All 
of the results reported in this paper constitute averages 
over tens of snapshots of the system (typically 50-100 
snapshots), with each snapshot being separated by an 
interval of the order of the eddy turnover time. 

C. Numerical iVIethod 

We solve equations (|ll2p and (|4l5p on a triply periodic 
domain using standard pseudospectral methods. The 
time-advancement of the diffusive terms is carried out 
exactly using the integrating factor method, while the 
remaining terms are treated using a third-order Runge- 
Kutta scheme. For a detailed description of the numerical 
method, sec, e.g., Ref. [3l|. We have conducted a number 
of MHD and RMHD simulations in both the balanced 
and imbalanced regimes. The parameters for each of the 
simulations are shown in Tabic HI 

III. RESULTS 

A. Comparing RMHD with MHD 

Here we directly compare the numerical solution of the 
full MHD system HUD with the RMHD system dUSD in 
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FIG. 1. A comparison of the MHD and RMHD field-perpendicular energy spectra (a) and the alignment angle (b) for balanced 
turbulence (Cases Ml and Rl). 



TABLE I. The parameters for each of the simulations. The 
first letter in the case number denotes the regime: M (MHD; 
eq. lU) or R (RMHD; eq. (|4])). In all cases v\\ — v, except for 
cases R2 and R3 where u\\ = 2.bu. 



Case 




iV|| 






Re = l/v 


(Tc 


Ml 


512 


512 


6 


1 


1800 





Rl 


512 


512 


6 


1 


1800 





M2 


1024 


256 


10 


1 


5600 


0.5 


R2 


1024 


256 


10 


1 


5600 


0.5 


M3 


256 


256 


6 


1 


800 





M4 


1024 


1024 


6 


1 


3200 





R3 


2048 


512 


10 


1 


14000 


0.5 


R4 


512 


256 


10 


1 


2200 


0.5 


M5 


512 


512 


5 


1 


1800 





M6 


512 


512 


5 


1 


3200 





M7 


512 


512 


5 


1 


4000 





M8 


1024 


1024 


5 


1 


4000 






both the balanced and imbalanced regimes. For the bal- 
anced case we compare simulations Ml and Rl. For the 
imbalanced case we compare the results of M2 and R2, 
for which the ratio of the Elsasser energies /E^ = 
(l + a,)/(l-a,)«3. 

Shown in Figure [T^ are the field-perpendicular en- 
ergy spectra for the velocity E'"{k±), the magnetic field 
E''{kj_) and the total spectrum E"^ = E"" + E'' for the 
case of balanced turbulence. Here 



(7) 



where q represents either the velocity or the magnetic 
field, q(fcj_) is the two-dimensional Fourier transfor- 
mation of q(x) in a plane perpendicular to Bq and 
k± = (fc^ -I- fcy)^/^. The average is taken over all field- 
perpendicular planes in the data cube and then over 
all data cubes (i.e. snapshots). The resulting field- 
perpendicular spectrum is equivalent to that obtained by 
integrating the three-dimensional Fourier spectrum over 
kz- It is clear that the three MHD and RMHD spectra 
are very similar beyond the forcing scales (the difference 



at k± = 1 , 2 is due to a slight difference in the forc- 
ing mechanisms in the MHD and RMHD cases, with the 
RMHD system ^ actually being solved for the Elsasser 
potentials 0^ , where = x V(/)^ , and hence incom- 
prcssibility being satisfied automatically, while for MHD 
the forces are constrained to be solenoidal). In particular, 
in both cases the total energy spectra are more closely 
matched with E'^(k±) oc kj^^^ than k^^^ , with the in- 
ertial range corresponding to 4 < fc_L < 20. It is also 
the case that for both MHD and RMHD the magnetic 
spectrum is slightly steeper than k^ and the velocity 
spectrum is slightly flatter. This interesting finding is 
pursued in more detail in Ref. |32| . 

Figure [TJd compares the alignment angle 6{r) between 
the shear Alfven velocity and magnetic fluctuations, 
which we deflne through the ratio 



%) 



((5Vr X (5br) 
(|<SVr||(5br|)' 



(8) 



Here (5vr = v(x) — v(x -I- r) and r is a point separa- 
tion vector in the plane perpendicular to the background 
magnetic field Bq. The average is taken over all points 
in a field-perpendicular slice, over all such slices in the 
data cube, and over all data cubes. In the MHD case the 
pseudo Alfven fluctuations are eliminated by removing 
the parts of (5vr and 5hj. that are in the direction of the 
local guide field, i.e. we construct 5vr = <^Vr — ((5vr • n)n 
where n = B(x)/|B(x)|. For RMHD the projection is 
not necessary. Figure [TJ) illustrates that the alignment 
angle in the MHD and RMHD regimes are almost in- 
distinguishable and that both are in excellent agreement 
with the theoretical prediction 6 oc r^^'^. 

A comparison of the ficld-pcrpcndicular Elsasser spec- 
tra £;^(fc_L) in the MHD and RMHD regimes of imbal- 
anced turbulence is shown in Figure [21 Again the agree- 
ment between the two systems is very good. In partic- 
ular, we find that weaker field obeys E~(k±) oc k_^_^^^, 
while E~^{k±) is slightly steeper at this resolution. In 
the next section we illustrate that the spectrum for the 
stronger field becomes flatter as the Reynolds number in- 
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FIG. 2. A comparison of the field-perpendicular energy spec- 
tra for the case of imbalanced MHD and RMHD turbu- 
lence (Cases M2 and R2). 



creases, implying that the universal regime has not yet 
been reached. Nonetheless, it is evident that the MHD 
and RMHD systems behave similarly (see also Ref. |33| ) . 

In summary, we find that in both the balanced and im- 
balanced regimes, the pseudo Alfven waves do not signif- 
icantly impact the strong turbulent cascade and RMHD 
accurately captures the MHD dynamics. 



B. Increasing the numerical resolution and Re 

We now investigate the robustness of the scalings laws 
as the numerical resolution and the Reynolds number 
increase and hence the extent of the incrtial range grows. 

Illustrated in Figure |3] are the results for a set of bal- 
anced MHD simulations corresponding to a four-fold in- 
crease in the numerical resolution. We note that all 
of the simulations in this set have a strong guide field, 
the computational domain is elongated in the z-direction 
in proportion to i?o, the number of mesh points in the 
field-parallel and field-perpendicular directions are equal, 
and the turbulence is excited in the strong state (i.e. 
the wavenumbers driven and the forcing correlation time 
are chosen so that the critical balance condition holds 
at the driving scales). Figure shows that the total 
field-perpendicular energy spectrum maintains the scal- 
ing E'^{k±) oc kj_^^'^ as the inertial range grows. We also 
note that the spectra fall off smoothly at large wavenum- 
bers, i.e. that the 'bottleneck effect' that leads to a 
pile up of energy at small scales is not present. Fig- 
ure [5)d illustrates that the scaling of the alignment angle 
is also in excellent agreement with the theoretical predic- 
tion 0r oc r^^^, with the extent of the region over which 
this scaling holds increasing as the resolution (and Re) 
increases. In fact, as explained in Ref. [l^ we believe 
that the alignment angle displays a significantly extended 
self-similar region that persists deep into the dissipation 
range, with the saturation of the r^^^ scaling being con- 
trolled by the smallest resolved scale of the simulation 



and hence decreasing by a factor of 2 as the resolution 
doubles. 

The effects on the Elsasser spectra E^{k±) of increas- 
ing the numerical resolution by a factor of four in the 
case of imbalanced RMHD turbulence are shown in Fig- 
ure [H It is seen that E^{k±) cx kj^^^ for all Reynolds 
numbers considered. However, the scaling of the stronger 
Elsasser field is difficult to pin down at these resolutions. 
Indeed, E^{k\^ appears to flatten as Re increases, and 
thus the universal regime has not yet been reached. In- 
deed, this is consistent with the Reynolds number for the 
stronger Elsasser field Re'^ ^ z~l/v being smaller than 
that for the weaker field. However, since the spectra are 
pinncd^i^ at the dissipation scale and anchored at the 
driving scale we proposed in Ref. [s^ that at sufficiently 
large Reynolds numbers the two spectra will become par- 
allel and obey the scaling E^{k±_) oc k^^"^ . 

C. Contamination of the scaling laws 

Having established above that for balanced turbulence 
the inertial range scalings are robust in optimally de- 
signed, carefully conducted numerical simulations, we 
now proceed to describe how numerical measurements of 
both the spectral exponent and the alignment angle can 
be contaminated. We emphasize that the scaling proper- 
ties are spoilt entirely due to numerical rather than phys- 
ical effects. The numerical errors can originate either at 
small scales and then back scatter to ultimately contam- 
inate the inertial range, or they can arise at large scales 
due to less than optimal numerical settings for studying 
strong magnetised turbulence. 

For strong balanced turbulence, it is known that nu- 
merical measurements of the total field-perpendicular 
energy spectrum can yield an exponent steeper than 
—3/2 when the properties of the driving are poorly con- 
trolled. For example, it was shown in Ref. |36| that, in 
the RMHD regime, as the number of modes driven in 
the field-parallel direction is increased the spectrum be- 
comes steeper, eventually approaching the weak turbu- 
lence result E{ki_) ~ fcjfjRefs. [13 andlH). A similar 
effect was found in Ref. [12 for MHD turbulence. Such 
effects are a result of the limited extent of the inertial 
range. Since present day computing power limits the in- 
ertial scales to approximately a decade in length, being 
able to reach the universal regime strongly relies on being 
able to limit the extent of the transition from the driving 
scales as much as possible. Such a contamination of the 
inertial scales would not arise if there was a significant 
separation between the driving scales and the small scale 
turbulence, as is the case in astrophysical settings where 
the Reynolds numbers are estimated to be orders of mag- 
nitude greater than the few thousands that present day 
computin g po wer permits. 

In Ref. Il6l . we explained how the alignment angle in 
balanced MHD turbulence can be spoilt by a strongly 
decreased field parallel resolution. When the simulation 
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FIG. 3. The total field-perpendicular energy spectrum (a) and the alignment angle (b) for a series of MHD calculations 
corresponding to increasing the numerical resolution and the Reynolds number (Cases M3, Ml and M4). 
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FIG. 4. The field-perpendicular Elsasser spectra for a series of 
RMHD calculations corresponding to increasing the numerical 
resolution and the Reynolds number (Cases R4, R2 and R3). 



FIG. 5. The alignment angle as the Reynolds number in- 
creases at fixed resolution (Cases M5, M6 and M7). 



is well resolved, alignment persists to much smaller scales 
than those over which the energy spectra display a power 
law behaviour. In Ref . [Tol we proposed that this could be 
due alignment being measured as the ratio of two struc- 
ture functions whose non- universal features cancel. Con- 
sequently, accurate measurement of the extended scaling 
behaviour relies on adequately resolving the small scale 
physics. Clearly this is spoilt when the resolution de- 
grades. Figure [5] illustrates that a similar contamina- 
tion occurs when the Reynolds number is pushed to the 
extreme. However, it is important that such behaviour 
should not be interpreted as a breakdown of the align- 
ment mechanism at high Re. Indeed, as is shown in Fig- 
ure [H excellent agreement with the r^/^ scaling can be 
recovered by increasing the numerical resolution. A simi- 
lar flattening of the alignment angle was shown in Ref. [2^ 
for a simulation with a reduced field-parallel resolution 
and sixth-order hyperviscosity. 




0.001 0.010 0.100 



r 

FIG. 6. The scaling of the alignment angle at high Re can be 
recovered by increasing the numerical resolution (Cases M7, 
M8). 



IV. DISCUSSION 

High-resolution direct numerical simulations play a key 
role in developing the theory for strong field-guided MHD 
turbulence. Success relies on working hard to harness all 
of the available computational power into reaching the 
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universal regime. When this is achieved, the numeri- 
cal data can be used to identify the preferred physical 
description when there exist conflicting theoretical mod- 
els. Furthermore, the simulations often reveal interesting 
new effects that subsequently guide further theoretical 
progress. 

We have discussed how the numerical setup can be 
optimised in order to accurately simulate strong MHD 
turbulence. By direct comparison with the full MHD 
system, we have also verified that the reduced MHD 
equations accurately model the turbulent cascade while 
being twice as computationally efficient. We have dis- 
cussed how balanced turbulence comprises local domains 
of highly aligned and anti-aligned velocity and magnetic 
fluctuations, obeying the theoretically predicted align- 
ment scaling 9 cx r^/"*. Consistent with the theory of 
dynamic alignment is the finding that the total field- 
perpendicular energy spectrum E{k±^) oc k^^^^ , with this 
scaling being maintained throughout a four-fold increase 
in numerical resolution. A similar set of simulations of 
imbalanced turbulence at steadily increasing Reynolds 
numbers up to Re = 14000 have shown that the weaker 
Elsasscr field obeys the same scaling, E~{kj_) cx kj_^^^. 
While a convincing result for the stronger Elsasscr field 
must await higher resolution tests, we propose that in 
the high Re limit the two spectra will become parallel 
and attain the scaling E^{kj_) cx kj^^^. 

We have also shown how sensitive numerical measure- 
ments can be to the design of the simulation. In particu- 
lar, we have shown how poorly resolved simulations will 
spoil measurements of the alignment angle and we have 
discussed how the spectral exponent can be sensitive to 
the physics of the large scale driving. The problems are 
due to the fact that computational power severely lim- 
its the extent of the inertial range, and as a consequence 
numerical effects can hamper our ability to infer the cor- 
rect physics. In the ideal case, any contamination at the 
driving or dissipative scales should not be felt sufficiently 
deep in the inertial range. Until enormous increases in 
computational power permit the investigation of such a 
regime, one must strive to ensure that the simulation is 
as clean as possible and that the numerical results are 
well converged. 
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